!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!

!--------------------------------------------------------------------------------------------------!
!   Copyright (c) 2008, 2009, Joost VandeVondele and Manuel Guidon                                 !
!   All rights reserved.                                                                           !
!                                                                                                  !
!   Redistribution and use in source and binary forms, with or without                             !
!   modification, are permitted provided that the following conditions are met:                    !
!       * Redistributions of source code must retain the above copyright                           !
!         notice, this list of conditions and the following disclaimer.                            !
!       * Redistributions in binary form must reproduce the above copyright                        !
!         notice, this list of conditions and the following disclaimer in the                      !
!         documentation and/or other materials provided with the distribution.                     !
!                                                                                                  !
!   THIS SOFTWARE IS PROVIDED BY Joost VandeVondele and Manuel Guidon AS IS AND ANY                !
!   EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED                      !
!   WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE                         !
!   DISCLAIMED. IN NO EVENT SHALL Joost VandeVondele or Manuel Guidon BE LIABLE FOR ANY            !
!   DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES                     !
!   (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;                   !
!   LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND                    !
!   ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT                     !
!   (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS                  !
!   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.                                   !
!--------------------------------------------------------------------------------------------------!

! **************************************************************************************************
!> \brief This module computes the basic integrals for the truncated coulomb operator
!>
!>          res(1) =G_0(R,T)= ((2*erf(sqrt(t))+erf(R-sqrt(t))-erf(R+sqrt(t)))/sqrt(t))
!>
!>        and up to 21 derivatives with respect to T
!>
!>          res(n+1)=(-1)**n d^n/dT^n G_0(R,T)
!>
!>        The function is only computed for values of R,T which fulfil
!>
!>          R**2 - 11.0_dp*R + 0.0_dp < T < R**2 + 11.0_dp*R + 50.0_dp where R>=0 T>=0
!>
!>        for T larger than the upper bound, 0 is returned
!>        (which is accurate at least up to 1.0E-16)
!>        while for T smaller than the lower bound, the caller is instructed
!>        to use the conventional gamma function instead
!>        (i.e. the limit of above expression for R to Infinity)
!>
!> \author Joost VandeVondele and Manuel Guidon
!> \par History
!>      Nov 2008, 2009 Joost VandeVondele and Manuel Guidon
!>      May 2019 A. Bussy: Added a get_maxl_init function to get current status of nderiv_init and
!>                moved the file to common (made it accessible from aobasis, same place as gamma.F).
! **************************************************************************************************
MODULE t_c_g0
   USE kinds,                           ONLY: dp
   USE message_passing,                 ONLY: mp_comm_type
#include "../base/base_uses.f90"

   IMPLICIT NONE
   PRIVATE

   PUBLIC :: t_c_g0_n, init, free_C0, get_lmax_init
   REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE, SAVE :: C0
   INTEGER, PARAMETER :: degree = 13
   REAL(KIND=dp), PARAMETER :: target_error = 0.100000E-08
   INTEGER, PARAMETER :: nderiv_max = 21
   INTEGER, SAVE      :: nderiv_init = -1
   INTEGER, SAVE      :: patches = -1

CONTAINS

! **************************************************************************************************
!> \brief ...
!> \param RES ...
!> \param use_gamma ...
!> \param R ...
!> \param T ...
!> \param NDERIV ...
! **************************************************************************************************
   SUBROUTINE t_c_g0_n(RES, use_gamma, R, T, NDERIV)
      REAL(KIND=dp), INTENT(OUT)                         :: RES(*)
      LOGICAL, INTENT(OUT)                               :: use_gamma
      REAL(KIND=dp), INTENT(IN)                          :: R, T
      INTEGER, INTENT(IN)                                :: NDERIV

      REAL(KIND=dp)                                      :: lower, TG1, TG2, upper, X1, X2

      use_gamma = .FALSE.
      upper = R**2 + 11.0_dp*R + 50.0_dp
      lower = R**2 - 11.0_dp*R + 0.0_dp
      IF (T > upper) THEN
         RES(1:NDERIV + 1) = 0.0_dp
         RETURN
      END IF
      IF (R <= 11.0_dp) THEN
         X2 = R/11.0_dp
         upper = R**2 + 11.0_dp*R + 50.0_dp
         lower = 0.0_dp
         X1 = (T - lower)/(upper - lower)
         IF (X1 <= 0.500000000000000000E+00_dp) THEN
            IF (X2 <= 0.500000000000000000E+00_dp) THEN
               IF (X2 <= 0.250000000000000000E+00_dp) THEN
                  IF (X2 <= 0.125000000000000000E+00_dp) THEN
                     IF (X1 <= 0.250000000000000000E+00_dp) THEN
                        IF (X2 <= 0.625000000000000000E-01_dp) THEN
                           IF (X1 <= 0.125000000000000000E+00_dp) THEN
                              IF (X2 <= 0.312500000000000000E-01_dp) THEN
                                 IF (X1 <= 0.625000000000000000E-01_dp) THEN
                                    IF (X2 <= 0.156250000000000000E-01_dp) THEN
                                       IF (X1 <= 0.312500000000000000E-01_dp) THEN
                                          TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                          TG2 = (2*X2 - 0.156250000000000000E-01_dp)*0.640000000000000000E+02_dp
                                          CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 1))
                                       ELSE
                                          TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                          TG2 = (2*X2 - 0.156250000000000000E-01_dp)*0.640000000000000000E+02_dp
                                          CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 2))
                                       END IF
                                    ELSE
                                       IF (X1 <= 0.312500000000000000E-01_dp) THEN
                                          TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                          TG2 = (2*X2 - 0.468750000000000000E-01_dp)*0.640000000000000000E+02_dp
                                          CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 3))
                                       ELSE
                                          TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                          TG2 = (2*X2 - 0.468750000000000000E-01_dp)*0.640000000000000000E+02_dp
                                          CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 4))
                                       END IF
                                    END IF
                                 ELSE
                                    IF (X2 <= 0.156250000000000000E-01_dp) THEN
                                       TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                       TG2 = (2*X2 - 0.156250000000000000E-01_dp)*0.640000000000000000E+02_dp
                                       CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 5))
                                    ELSE
                                       TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                       TG2 = (2*X2 - 0.468750000000000000E-01_dp)*0.640000000000000000E+02_dp
                                       CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 6))
                                    END IF
                                 END IF
                              ELSE
                                 IF (X1 <= 0.625000000000000000E-01_dp) THEN
                                    IF (X2 <= 0.468750000000000000E-01_dp) THEN
                                       IF (X1 <= 0.312500000000000000E-01_dp) THEN
                                          TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                          TG2 = (2*X2 - 0.781250000000000000E-01_dp)*0.640000000000000000E+02_dp
                                          CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 7))
                                       ELSE
                                          TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                          TG2 = (2*X2 - 0.781250000000000000E-01_dp)*0.640000000000000000E+02_dp
                                          CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 8))
                                       END IF
                                    ELSE
                                       IF (X1 <= 0.312500000000000000E-01_dp) THEN
                                          TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                          TG2 = (2*X2 - 0.109375000000000000E+00_dp)*0.640000000000000000E+02_dp
                                          CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 9))
                                       ELSE
                                          TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                          TG2 = (2*X2 - 0.109375000000000000E+00_dp)*0.640000000000000000E+02_dp
                                          CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 10))
                                       END IF
                                    END IF
                                 ELSE
                                    IF (X2 <= 0.468750000000000000E-01_dp) THEN
                                       TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                       TG2 = (2*X2 - 0.781250000000000000E-01_dp)*0.640000000000000000E+02_dp
                                       CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 11))
                                    ELSE
                                       TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                       TG2 = (2*X2 - 0.109375000000000000E+00_dp)*0.640000000000000000E+02_dp
                                       CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 12))
                                    END IF
                                 END IF
                              END IF
                           ELSE
                              IF (X2 <= 0.312500000000000000E-01_dp) THEN
                                 IF (X1 <= 0.187500000000000000E+00_dp) THEN
                                    TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                    TG2 = (2*X2 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                    CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 13))
                                 ELSE
                                    TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                    TG2 = (2*X2 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                    CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 14))
                                 END IF
                              ELSE
                                 IF (X1 <= 0.187500000000000000E+00_dp) THEN
                                    TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                    TG2 = (2*X2 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                    CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 15))
                                 ELSE
                                    TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                    TG2 = (2*X2 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                    CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 16))
                                 END IF
                              END IF
                           END IF
                        ELSE
                           IF (X1 <= 0.125000000000000000E+00_dp) THEN
                              IF (X2 <= 0.937500000000000000E-01_dp) THEN
                                 IF (X1 <= 0.625000000000000000E-01_dp) THEN
                                    IF (X2 <= 0.781250000000000000E-01_dp) THEN
                                       IF (X1 <= 0.312500000000000000E-01_dp) THEN
                                          TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                          TG2 = (2*X2 - 0.140625000000000000E+00_dp)*0.640000000000000000E+02_dp
                                          CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 17))
                                       ELSE
                                          TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                          TG2 = (2*X2 - 0.140625000000000000E+00_dp)*0.640000000000000000E+02_dp
                                          CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 18))
                                       END IF
                                    ELSE
                                       IF (X1 <= 0.312500000000000000E-01_dp) THEN
                                          TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                          TG2 = (2*X2 - 0.171875000000000000E+00_dp)*0.640000000000000000E+02_dp
                                          CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 19))
                                       ELSE
                                          TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                          TG2 = (2*X2 - 0.171875000000000000E+00_dp)*0.640000000000000000E+02_dp
                                          CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 20))
                                       END IF
                                    END IF
                                 ELSE
                                    IF (X2 <= 0.781250000000000000E-01_dp) THEN
                                       TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                       TG2 = (2*X2 - 0.140625000000000000E+00_dp)*0.640000000000000000E+02_dp
                                       CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 21))
                                    ELSE
                                       TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                       TG2 = (2*X2 - 0.171875000000000000E+00_dp)*0.640000000000000000E+02_dp
                                       CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 22))
                                    END IF
                                 END IF
                              ELSE
                                 IF (X1 <= 0.625000000000000000E-01_dp) THEN
                                    IF (X2 <= 0.109375000000000000E+00_dp) THEN
                                       IF (X1 <= 0.312500000000000000E-01_dp) THEN
                                          TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                          TG2 = (2*X2 - 0.203125000000000000E+00_dp)*0.640000000000000000E+02_dp
                                          CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 23))
                                       ELSE
                                          TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                          TG2 = (2*X2 - 0.203125000000000000E+00_dp)*0.640000000000000000E+02_dp
                                          CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 24))
                                       END IF
                                    ELSE
                                       IF (X1 <= 0.312500000000000000E-01_dp) THEN
                                          TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                          TG2 = (2*X2 - 0.234375000000000000E+00_dp)*0.640000000000000000E+02_dp
                                          CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 25))
                                       ELSE
                                          TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                          TG2 = (2*X2 - 0.234375000000000000E+00_dp)*0.640000000000000000E+02_dp
                                          CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 26))
                                       END IF
                                    END IF
                                 ELSE
                                    IF (X2 <= 0.109375000000000000E+00_dp) THEN
                                       TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                       TG2 = (2*X2 - 0.203125000000000000E+00_dp)*0.640000000000000000E+02_dp
                                       CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 27))
                                    ELSE
                                       TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                       TG2 = (2*X2 - 0.234375000000000000E+00_dp)*0.640000000000000000E+02_dp
                                       CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 28))
                                    END IF
                                 END IF
                              END IF
                           ELSE
                              IF (X1 <= 0.187500000000000000E+00_dp) THEN
                                 IF (X2 <= 0.937500000000000000E-01_dp) THEN
                                    TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                    TG2 = (2*X2 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
                                    CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 29))
                                 ELSE
                                    TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                    TG2 = (2*X2 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
                                    CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 30))
                                 END IF
                              ELSE
                                 TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 31))
                              END IF
                           END IF
                        END IF
                     ELSE
                        IF (X1 <= 0.375000000000000000E+00_dp) THEN
                           TG1 = (2*X1 - 0.625000000000000000E+00_dp)*0.800000000000000000E+01_dp
                           TG2 = (2*X2 - 0.125000000000000000E+00_dp)*0.800000000000000000E+01_dp
                           CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 32))
                        ELSE
                           TG1 = (2*X1 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
                           TG2 = (2*X2 - 0.125000000000000000E+00_dp)*0.800000000000000000E+01_dp
                           CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 33))
                        END IF
                     END IF
                  ELSE
                     IF (X1 <= 0.250000000000000000E+00_dp) THEN
                        IF (X2 <= 0.187500000000000000E+00_dp) THEN
                           IF (X1 <= 0.125000000000000000E+00_dp) THEN
                              IF (X2 <= 0.156250000000000000E+00_dp) THEN
                                 IF (X1 <= 0.625000000000000000E-01_dp) THEN
                                    IF (X1 <= 0.312500000000000000E-01_dp) THEN
                                       IF (X2 <= 0.140625000000000000E+00_dp) THEN
                                          TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                          TG2 = (2*X2 - 0.265625000000000000E+00_dp)*0.640000000000000000E+02_dp
                                          CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 34))
                                       ELSE
                                          TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                          TG2 = (2*X2 - 0.296875000000000000E+00_dp)*0.640000000000000000E+02_dp
                                          CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 35))
                                       END IF
                                    ELSE
                                       TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                       TG2 = (2*X2 - 0.281250000000000000E+00_dp)*0.320000000000000000E+02_dp
                                       CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 36))
                                    END IF
                                 ELSE
                                    IF (X1 <= 0.937500000000000000E-01_dp) THEN
                                       TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
                                       TG2 = (2*X2 - 0.281250000000000000E+00_dp)*0.320000000000000000E+02_dp
                                       CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 37))
                                    ELSE
                                       TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
                                       TG2 = (2*X2 - 0.281250000000000000E+00_dp)*0.320000000000000000E+02_dp
                                       CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 38))
                                    END IF
                                 END IF
                              ELSE
                                 IF (X1 <= 0.625000000000000000E-01_dp) THEN
                                    IF (X1 <= 0.312500000000000000E-01_dp) THEN
                                       IF (X2 <= 0.171875000000000000E+00_dp) THEN
                                          TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                          TG2 = (2*X2 - 0.328125000000000000E+00_dp)*0.640000000000000000E+02_dp
                                          CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 39))
                                       ELSE
                                          TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                          TG2 = (2*X2 - 0.359375000000000000E+00_dp)*0.640000000000000000E+02_dp
                                          CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 40))
                                       END IF
                                    ELSE
                                       TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                       TG2 = (2*X2 - 0.343750000000000000E+00_dp)*0.320000000000000000E+02_dp
                                       CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 41))
                                    END IF
                                 ELSE
                                    TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                    TG2 = (2*X2 - 0.343750000000000000E+00_dp)*0.320000000000000000E+02_dp
                                    CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 42))
                                 END IF
                              END IF
                           ELSE
                              IF (X1 <= 0.187500000000000000E+00_dp) THEN
                                 TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 43))
                              ELSE
                                 TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 44))
                              END IF
                           END IF
                        ELSE
                           IF (X1 <= 0.125000000000000000E+00_dp) THEN
                              IF (X1 <= 0.625000000000000000E-01_dp) THEN
                                 IF (X2 <= 0.218750000000000000E+00_dp) THEN
                                    IF (X1 <= 0.312500000000000000E-01_dp) THEN
                                       IF (X2 <= 0.203125000000000000E+00_dp) THEN
                                          TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                          TG2 = (2*X2 - 0.390625000000000000E+00_dp)*0.640000000000000000E+02_dp
                                          CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 45))
                                       ELSE
                                          TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                          TG2 = (2*X2 - 0.421875000000000000E+00_dp)*0.640000000000000000E+02_dp
                                          CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 46))
                                       END IF
                                    ELSE
                                       TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                       TG2 = (2*X2 - 0.406250000000000000E+00_dp)*0.320000000000000000E+02_dp
                                       CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 47))
                                    END IF
                                 ELSE
                                    IF (X1 <= 0.312500000000000000E-01_dp) THEN
                                       IF (X2 <= 0.234375000000000000E+00_dp) THEN
                                          TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                          TG2 = (2*X2 - 0.453125000000000000E+00_dp)*0.640000000000000000E+02_dp
                                          CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 48))
                                       ELSE
                                          TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                          TG2 = (2*X2 - 0.484375000000000000E+00_dp)*0.640000000000000000E+02_dp
                                          CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 49))
                                       END IF
                                    ELSE
                                       TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                       TG2 = (2*X2 - 0.468750000000000000E+00_dp)*0.320000000000000000E+02_dp
                                       CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 50))
                                    END IF
                                 END IF
                              ELSE
                                 IF (X2 <= 0.218750000000000000E+00_dp) THEN
                                    TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                    TG2 = (2*X2 - 0.406250000000000000E+00_dp)*0.320000000000000000E+02_dp
                                    CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 51))
                                 ELSE
                                    TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                    TG2 = (2*X2 - 0.468750000000000000E+00_dp)*0.320000000000000000E+02_dp
                                    CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 52))
                                 END IF
                              END IF
                           ELSE
                              IF (X1 <= 0.187500000000000000E+00_dp) THEN
                                 TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 53))
                              ELSE
                                 TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 54))
                              END IF
                           END IF
                        END IF
                     ELSE
                        IF (X1 <= 0.375000000000000000E+00_dp) THEN
                           IF (X1 <= 0.312500000000000000E+00_dp) THEN
                              TG1 = (2*X1 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
                              TG2 = (2*X2 - 0.375000000000000000E+00_dp)*0.800000000000000000E+01_dp
                              CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 55))
                           ELSE
                              TG1 = (2*X1 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
                              TG2 = (2*X2 - 0.375000000000000000E+00_dp)*0.800000000000000000E+01_dp
                              CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 56))
                           END IF
                        ELSE
                           TG1 = (2*X1 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
                           TG2 = (2*X2 - 0.375000000000000000E+00_dp)*0.800000000000000000E+01_dp
                           CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 57))
                        END IF
                     END IF
                  END IF
               ELSE
                  IF (X1 <= 0.250000000000000000E+00_dp) THEN
                     IF (X1 <= 0.125000000000000000E+00_dp) THEN
                        IF (X1 <= 0.625000000000000000E-01_dp) THEN
                           IF (X2 <= 0.375000000000000000E+00_dp) THEN
                              IF (X2 <= 0.312500000000000000E+00_dp) THEN
                                 IF (X1 <= 0.312500000000000000E-01_dp) THEN
                                    IF (X2 <= 0.281250000000000000E+00_dp) THEN
                                       TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                       TG2 = (2*X2 - 0.531250000000000000E+00_dp)*0.320000000000000000E+02_dp
                                       CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 58))
                                    ELSE
                                       TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                       TG2 = (2*X2 - 0.593750000000000000E+00_dp)*0.320000000000000000E+02_dp
                                       CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 59))
                                    END IF
                                 ELSE
                                    IF (X2 <= 0.281250000000000000E+00_dp) THEN
                                       TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                       TG2 = (2*X2 - 0.531250000000000000E+00_dp)*0.320000000000000000E+02_dp
                                       CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 60))
                                    ELSE
                                       TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                       TG2 = (2*X2 - 0.593750000000000000E+00_dp)*0.320000000000000000E+02_dp
                                       CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 61))
                                    END IF
                                 END IF
                              ELSE
                                 IF (X1 <= 0.312500000000000000E-01_dp) THEN
                                    IF (X2 <= 0.343750000000000000E+00_dp) THEN
                                       TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                       TG2 = (2*X2 - 0.656250000000000000E+00_dp)*0.320000000000000000E+02_dp
                                       CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 62))
                                    ELSE
                                       TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                       TG2 = (2*X2 - 0.718750000000000000E+00_dp)*0.320000000000000000E+02_dp
                                       CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 63))
                                    END IF
                                 ELSE
                                    TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                    TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                    CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 64))
                                 END IF
                              END IF
                           ELSE
                              IF (X1 <= 0.312500000000000000E-01_dp) THEN
                                 IF (X2 <= 0.437500000000000000E+00_dp) THEN
                                    IF (X1 <= 0.156250000000000000E-01_dp) THEN
                                       TG1 = (2*X1 - 0.156250000000000000E-01_dp)*0.640000000000000000E+02_dp
                                       TG2 = (2*X2 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                       CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 65))
                                    ELSE
                                       TG1 = (2*X1 - 0.468750000000000000E-01_dp)*0.640000000000000000E+02_dp
                                       TG2 = (2*X2 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                       CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 66))
                                    END IF
                                 ELSE
                                    IF (X1 <= 0.156250000000000000E-01_dp) THEN
                                       TG1 = (2*X1 - 0.156250000000000000E-01_dp)*0.640000000000000000E+02_dp
                                       TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                       CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 67))
                                    ELSE
                                       TG1 = (2*X1 - 0.468750000000000000E-01_dp)*0.640000000000000000E+02_dp
                                       TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                       CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 68))
                                    END IF
                                 END IF
                              ELSE
                                 IF (X2 <= 0.437500000000000000E+00_dp) THEN
                                    TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                    TG2 = (2*X2 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                    CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 69))
                                 ELSE
                                    TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                    TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                    CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 70))
                                 END IF
                              END IF
                           END IF
                        ELSE
                           IF (X2 <= 0.375000000000000000E+00_dp) THEN
                              IF (X2 <= 0.312500000000000000E+00_dp) THEN
                                 IF (X1 <= 0.937500000000000000E-01_dp) THEN
                                    TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
                                    TG2 = (2*X2 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                    CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 71))
                                 ELSE
                                    TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
                                    TG2 = (2*X2 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                    CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 72))
                                 END IF
                              ELSE
                                 IF (X1 <= 0.937500000000000000E-01_dp) THEN
                                    TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
                                    TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                    CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 73))
                                 ELSE
                                    TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
                                    TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                    CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 74))
                                 END IF
                              END IF
                           ELSE
                              IF (X1 <= 0.937500000000000000E-01_dp) THEN
                                 IF (X2 <= 0.437500000000000000E+00_dp) THEN
                                    TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
                                    TG2 = (2*X2 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                    CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 75))
                                 ELSE
                                    TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
                                    TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                    CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 76))
                                 END IF
                              ELSE
                                 IF (X2 <= 0.437500000000000000E+00_dp) THEN
                                    TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
                                    TG2 = (2*X2 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                    CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 77))
                                 ELSE
                                    TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
                                    TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                    CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 78))
                                 END IF
                              END IF
                           END IF
                        END IF
                     ELSE
                        IF (X2 <= 0.375000000000000000E+00_dp) THEN
                           IF (X1 <= 0.187500000000000000E+00_dp) THEN
                              IF (X2 <= 0.312500000000000000E+00_dp) THEN
                                 TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 79))
                              ELSE
                                 IF (X1 <= 0.156250000000000000E+00_dp) THEN
                                    TG1 = (2*X1 - 0.281250000000000000E+00_dp)*0.320000000000000000E+02_dp
                                    TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                    CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 80))
                                 ELSE
                                    TG1 = (2*X1 - 0.343750000000000000E+00_dp)*0.320000000000000000E+02_dp
                                    TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                    CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 81))
                                 END IF
                              END IF
                           ELSE
                              IF (X2 <= 0.312500000000000000E+00_dp) THEN
                                 TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 82))
                              ELSE
                                 TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 83))
                              END IF
                           END IF
                        ELSE
                           IF (X1 <= 0.187500000000000000E+00_dp) THEN
                              IF (X2 <= 0.437500000000000000E+00_dp) THEN
                                 TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 84))
                              ELSE
                                 IF (X1 <= 0.156250000000000000E+00_dp) THEN
                                    TG1 = (2*X1 - 0.281250000000000000E+00_dp)*0.320000000000000000E+02_dp
                                    TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                    CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 85))
                                 ELSE
                                    TG1 = (2*X1 - 0.343750000000000000E+00_dp)*0.320000000000000000E+02_dp
                                    TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
                                    CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 86))
                                 END IF
                              END IF
                           ELSE
                              IF (X1 <= 0.218750000000000000E+00_dp) THEN
                                 TG1 = (2*X1 - 0.406250000000000000E+00_dp)*0.320000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 87))
                              ELSE
                                 TG1 = (2*X1 - 0.468750000000000000E+00_dp)*0.320000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 88))
                              END IF
                           END IF
                        END IF
                     END IF
                  ELSE
                     IF (X1 <= 0.375000000000000000E+00_dp) THEN
                        IF (X2 <= 0.375000000000000000E+00_dp) THEN
                           IF (X1 <= 0.312500000000000000E+00_dp) THEN
                              TG1 = (2*X1 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
                              TG2 = (2*X2 - 0.625000000000000000E+00_dp)*0.800000000000000000E+01_dp
                              CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 89))
                           ELSE
                              TG1 = (2*X1 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
                              TG2 = (2*X2 - 0.625000000000000000E+00_dp)*0.800000000000000000E+01_dp
                              CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 90))
                           END IF
                        ELSE
                           IF (X1 <= 0.312500000000000000E+00_dp) THEN
                              IF (X1 <= 0.281250000000000000E+00_dp) THEN
                                 TG1 = (2*X1 - 0.531250000000000000E+00_dp)*0.320000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 91))
                              ELSE
                                 TG1 = (2*X1 - 0.593750000000000000E+00_dp)*0.320000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 92))
                              END IF
                           ELSE
                              TG1 = (2*X1 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
                              TG2 = (2*X2 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
                              CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 93))
                           END IF
                        END IF
                     ELSE
                        IF (X1 <= 0.437500000000000000E+00_dp) THEN
                           TG1 = (2*X1 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
                           TG2 = (2*X2 - 0.750000000000000000E+00_dp)*0.400000000000000000E+01_dp
                           CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 94))
                        ELSE
                           TG1 = (2*X1 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
                           TG2 = (2*X2 - 0.750000000000000000E+00_dp)*0.400000000000000000E+01_dp
                           CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 95))
                        END IF
                     END IF
                  END IF
               END IF
            ELSE
               IF (X1 <= 0.250000000000000000E+00_dp) THEN
                  IF (X1 <= 0.125000000000000000E+00_dp) THEN
                     IF (X1 <= 0.625000000000000000E-01_dp) THEN
                        IF (X1 <= 0.312500000000000000E-01_dp) THEN
                           IF (X1 <= 0.156250000000000000E-01_dp) THEN
                              IF (X1 <= 0.781250000000000000E-02_dp) THEN
                                 IF (X2 <= 0.750000000000000000E+00_dp) THEN
                                    TG1 = (2*X1 - 0.781250000000000000E-02_dp)*0.128000000000000000E+03_dp
                                    TG2 = (2*X2 - 0.125000000000000000E+01_dp)*0.400000000000000000E+01_dp
                                    CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 96))
                                 ELSE
                                    TG1 = (2*X1 - 0.781250000000000000E-02_dp)*0.128000000000000000E+03_dp
                                    TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
                                    CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 97))
                                 END IF
                              ELSE
                                 IF (X2 <= 0.750000000000000000E+00_dp) THEN
                                    TG1 = (2*X1 - 0.234375000000000000E-01_dp)*0.128000000000000000E+03_dp
                                    TG2 = (2*X2 - 0.125000000000000000E+01_dp)*0.400000000000000000E+01_dp
                                    CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 98))
                                 ELSE
                                    TG1 = (2*X1 - 0.234375000000000000E-01_dp)*0.128000000000000000E+03_dp
                                    TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
                                    CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 99))
                                 END IF
                              END IF
                           ELSE
                              IF (X2 <= 0.750000000000000000E+00_dp) THEN
                                 TG1 = (2*X1 - 0.468750000000000000E-01_dp)*0.640000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.125000000000000000E+01_dp)*0.400000000000000000E+01_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 100))
                              ELSE
                                 TG1 = (2*X1 - 0.468750000000000000E-01_dp)*0.640000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 101))
                              END IF
                           END IF
                        ELSE
                           IF (X2 <= 0.750000000000000000E+00_dp) THEN
                              IF (X2 <= 0.625000000000000000E+00_dp) THEN
                                 TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 102))
                              ELSE
                                 TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 103))
                              END IF
                           ELSE
                              IF (X1 <= 0.468750000000000000E-01_dp) THEN
                                 TG1 = (2*X1 - 0.781250000000000000E-01_dp)*0.640000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 104))
                              ELSE
                                 TG1 = (2*X1 - 0.109375000000000000E+00_dp)*0.640000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 105))
                              END IF
                           END IF
                        END IF
                     ELSE
                        IF (X2 <= 0.750000000000000000E+00_dp) THEN
                           IF (X2 <= 0.625000000000000000E+00_dp) THEN
                              IF (X1 <= 0.937500000000000000E-01_dp) THEN
                                 TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 106))
                              ELSE
                                 TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 107))
                              END IF
                           ELSE
                              IF (X1 <= 0.937500000000000000E-01_dp) THEN
                                 TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 108))
                              ELSE
                                 TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 109))
                              END IF
                           END IF
                        ELSE
                           IF (X1 <= 0.937500000000000000E-01_dp) THEN
                              TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
                              TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
                              CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 110))
                           ELSE
                              TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
                              TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
                              CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 111))
                           END IF
                        END IF
                     END IF
                  ELSE
                     IF (X2 <= 0.750000000000000000E+00_dp) THEN
                        IF (X2 <= 0.625000000000000000E+00_dp) THEN
                           IF (X1 <= 0.187500000000000000E+00_dp) THEN
                              IF (X1 <= 0.156250000000000000E+00_dp) THEN
                                 TG1 = (2*X1 - 0.281250000000000000E+00_dp)*0.320000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 112))
                              ELSE
                                 TG1 = (2*X1 - 0.343750000000000000E+00_dp)*0.320000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 113))
                              END IF
                           ELSE
                              IF (X1 <= 0.218750000000000000E+00_dp) THEN
                                 TG1 = (2*X1 - 0.406250000000000000E+00_dp)*0.320000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 114))
                              ELSE
                                 TG1 = (2*X1 - 0.468750000000000000E+00_dp)*0.320000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 115))
                              END IF
                           END IF
                        ELSE
                           IF (X1 <= 0.187500000000000000E+00_dp) THEN
                              IF (X1 <= 0.156250000000000000E+00_dp) THEN
                                 TG1 = (2*X1 - 0.281250000000000000E+00_dp)*0.320000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 116))
                              ELSE
                                 TG1 = (2*X1 - 0.343750000000000000E+00_dp)*0.320000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 117))
                              END IF
                           ELSE
                              IF (X1 <= 0.218750000000000000E+00_dp) THEN
                                 TG1 = (2*X1 - 0.406250000000000000E+00_dp)*0.320000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 118))
                              ELSE
                                 TG1 = (2*X1 - 0.468750000000000000E+00_dp)*0.320000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 119))
                              END IF
                           END IF
                        END IF
                     ELSE
                        IF (X1 <= 0.187500000000000000E+00_dp) THEN
                           IF (X2 <= 0.875000000000000000E+00_dp) THEN
                              TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
                              TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
                              CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 120))
                           ELSE
                              TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
                              TG2 = (2*X2 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
                              CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 121))
                           END IF
                        ELSE
                           IF (X2 <= 0.875000000000000000E+00_dp) THEN
                              IF (X1 <= 0.218750000000000000E+00_dp) THEN
                                 TG1 = (2*X1 - 0.406250000000000000E+00_dp)*0.320000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 122))
                              ELSE
                                 TG1 = (2*X1 - 0.468750000000000000E+00_dp)*0.320000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 123))
                              END IF
                           ELSE
                              TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
                              TG2 = (2*X2 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
                              CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 124))
                           END IF
                        END IF
                     END IF
                  END IF
               ELSE
                  IF (X1 <= 0.375000000000000000E+00_dp) THEN
                     IF (X2 <= 0.750000000000000000E+00_dp) THEN
                        IF (X1 <= 0.312500000000000000E+00_dp) THEN
                           IF (X2 <= 0.625000000000000000E+00_dp) THEN
                              IF (X1 <= 0.281250000000000000E+00_dp) THEN
                                 TG1 = (2*X1 - 0.531250000000000000E+00_dp)*0.320000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 125))
                              ELSE
                                 TG1 = (2*X1 - 0.593750000000000000E+00_dp)*0.320000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 126))
                              END IF
                           ELSE
                              IF (X1 <= 0.281250000000000000E+00_dp) THEN
                                 TG1 = (2*X1 - 0.531250000000000000E+00_dp)*0.320000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 127))
                              ELSE
                                 TG1 = (2*X1 - 0.593750000000000000E+00_dp)*0.320000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 128))
                              END IF
                           END IF
                        ELSE
                           IF (X2 <= 0.625000000000000000E+00_dp) THEN
                              TG1 = (2*X1 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
                              TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
                              CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 129))
                           ELSE
                              IF (X1 <= 0.343750000000000000E+00_dp) THEN
                                 TG1 = (2*X1 - 0.656250000000000000E+00_dp)*0.320000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 130))
                              ELSE
                                 TG1 = (2*X1 - 0.718750000000000000E+00_dp)*0.320000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 131))
                              END IF
                           END IF
                        END IF
                     ELSE
                        IF (X1 <= 0.312500000000000000E+00_dp) THEN
                           IF (X2 <= 0.875000000000000000E+00_dp) THEN
                              IF (X1 <= 0.281250000000000000E+00_dp) THEN
                                 TG1 = (2*X1 - 0.531250000000000000E+00_dp)*0.320000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 132))
                              ELSE
                                 TG1 = (2*X1 - 0.593750000000000000E+00_dp)*0.320000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 133))
                              END IF
                           ELSE
                              IF (X1 <= 0.281250000000000000E+00_dp) THEN
                                 TG1 = (2*X1 - 0.531250000000000000E+00_dp)*0.320000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 134))
                              ELSE
                                 TG1 = (2*X1 - 0.593750000000000000E+00_dp)*0.320000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 135))
                              END IF
                           END IF
                        ELSE
                           IF (X2 <= 0.875000000000000000E+00_dp) THEN
                              IF (X1 <= 0.343750000000000000E+00_dp) THEN
                                 TG1 = (2*X1 - 0.656250000000000000E+00_dp)*0.320000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 136))
                              ELSE
                                 TG1 = (2*X1 - 0.718750000000000000E+00_dp)*0.320000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 137))
                              END IF
                           ELSE
                              TG1 = (2*X1 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
                              TG2 = (2*X2 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
                              CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 138))
                           END IF
                        END IF
                     END IF
                  ELSE
                     IF (X2 <= 0.750000000000000000E+00_dp) THEN
                        IF (X1 <= 0.437500000000000000E+00_dp) THEN
                           IF (X2 <= 0.625000000000000000E+00_dp) THEN
                              TG1 = (2*X1 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
                              TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
                              CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 139))
                           ELSE
                              TG1 = (2*X1 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
                              TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
                              CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 140))
                           END IF
                        ELSE
                           TG1 = (2*X1 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
                           TG2 = (2*X2 - 0.125000000000000000E+01_dp)*0.400000000000000000E+01_dp
                           CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 141))
                        END IF
                     ELSE
                        IF (X1 <= 0.437500000000000000E+00_dp) THEN
                           IF (X2 <= 0.875000000000000000E+00_dp) THEN
                              TG1 = (2*X1 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
                              TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
                              CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 142))
                           ELSE
                              TG1 = (2*X1 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
                              TG2 = (2*X2 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
                              CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 143))
                           END IF
                        ELSE
                           IF (X2 <= 0.875000000000000000E+00_dp) THEN
                              TG1 = (2*X1 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
                              TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
                              CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 144))
                           ELSE
                              TG1 = (2*X1 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
                              TG2 = (2*X2 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
                              CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 145))
                           END IF
                        END IF
                     END IF
                  END IF
               END IF
            END IF
         ELSE
            IF (X1 <= 0.750000000000000000E+00_dp) THEN
               IF (X2 <= 0.500000000000000000E+00_dp) THEN
                  IF (X1 <= 0.625000000000000000E+00_dp) THEN
                     IF (X2 <= 0.250000000000000000E+00_dp) THEN
                        TG1 = (2*X1 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
                        TG2 = (2*X2 - 0.250000000000000000E+00_dp)*0.400000000000000000E+01_dp
                        CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 146))
                     ELSE
                        TG1 = (2*X1 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
                        TG2 = (2*X2 - 0.750000000000000000E+00_dp)*0.400000000000000000E+01_dp
                        CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 147))
                     END IF
                  ELSE
                     TG1 = (2*X1 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
                     TG2 = (2*X2 - 0.500000000000000000E+00_dp)*0.200000000000000000E+01_dp
                     CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 148))
                  END IF
               ELSE
                  IF (X1 <= 0.625000000000000000E+00_dp) THEN
                     IF (X2 <= 0.750000000000000000E+00_dp) THEN
                        IF (X1 <= 0.562500000000000000E+00_dp) THEN
                           TG1 = (2*X1 - 0.106250000000000000E+01_dp)*0.160000000000000000E+02_dp
                           TG2 = (2*X2 - 0.125000000000000000E+01_dp)*0.400000000000000000E+01_dp
                           CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 149))
                        ELSE
                           TG1 = (2*X1 - 0.118750000000000000E+01_dp)*0.160000000000000000E+02_dp
                           TG2 = (2*X2 - 0.125000000000000000E+01_dp)*0.400000000000000000E+01_dp
                           CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 150))
                        END IF
                     ELSE
                        IF (X1 <= 0.562500000000000000E+00_dp) THEN
                           TG1 = (2*X1 - 0.106250000000000000E+01_dp)*0.160000000000000000E+02_dp
                           TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
                           CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 151))
                        ELSE
                           TG1 = (2*X1 - 0.118750000000000000E+01_dp)*0.160000000000000000E+02_dp
                           TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
                           CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 152))
                        END IF
                     END IF
                  ELSE
                     IF (X1 <= 0.687500000000000000E+00_dp) THEN
                        TG1 = (2*X1 - 0.131250000000000000E+01_dp)*0.160000000000000000E+02_dp
                        TG2 = (2*X2 - 0.150000000000000000E+01_dp)*0.200000000000000000E+01_dp
                        CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 153))
                     ELSE
                        TG1 = (2*X1 - 0.143750000000000000E+01_dp)*0.160000000000000000E+02_dp
                        TG2 = (2*X2 - 0.150000000000000000E+01_dp)*0.200000000000000000E+01_dp
                        CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 154))
                     END IF
                  END IF
               END IF
            ELSE
               TG1 = (2*X1 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
               TG2 = (2*X2 - 0.100000000000000000E+01_dp)*0.100000000000000000E+01_dp
               CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 155))
            END IF
         END IF
      ELSE
         IF (T < lower) THEN
            use_gamma = .TRUE.
            RETURN
         END IF
         X2 = 11.0_dp/R
         X1 = (T - lower)/(upper - lower)
         IF (X1 <= 0.500000000000000000E+00_dp) THEN
            IF (X1 <= 0.250000000000000000E+00_dp) THEN
               IF (X2 <= 0.500000000000000000E+00_dp) THEN
                  IF (X1 <= 0.125000000000000000E+00_dp) THEN
                     IF (X2 <= 0.250000000000000000E+00_dp) THEN
                        TG1 = (2*X1 - 0.125000000000000000E+00_dp)*0.800000000000000000E+01_dp
                        TG2 = (2*X2 - 0.250000000000000000E+00_dp)*0.400000000000000000E+01_dp
                        CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 156))
                     ELSE
                        TG1 = (2*X1 - 0.125000000000000000E+00_dp)*0.800000000000000000E+01_dp
                        TG2 = (2*X2 - 0.750000000000000000E+00_dp)*0.400000000000000000E+01_dp
                        CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 157))
                     END IF
                  ELSE
                     TG1 = (2*X1 - 0.375000000000000000E+00_dp)*0.800000000000000000E+01_dp
                     TG2 = (2*X2 - 0.500000000000000000E+00_dp)*0.200000000000000000E+01_dp
                     CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 158))
                  END IF
               ELSE
                  IF (X1 <= 0.125000000000000000E+00_dp) THEN
                     IF (X2 <= 0.750000000000000000E+00_dp) THEN
                        IF (X2 <= 0.625000000000000000E+00_dp) THEN
                           TG1 = (2*X1 - 0.125000000000000000E+00_dp)*0.800000000000000000E+01_dp
                           TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
                           CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 159))
                        ELSE
                           IF (X1 <= 0.625000000000000000E-01_dp) THEN
                              TG1 = (2*X1 - 0.625000000000000000E-01_dp)*0.160000000000000000E+02_dp
                              TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
                              CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 160))
                           ELSE
                              TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
                              TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
                              CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 161))
                           END IF
                        END IF
                     ELSE
                        IF (X1 <= 0.625000000000000000E-01_dp) THEN
                           IF (X2 <= 0.875000000000000000E+00_dp) THEN
                              IF (X1 <= 0.312500000000000000E-01_dp) THEN
                                 IF (X2 <= 0.812500000000000000E+00_dp) THEN
                                    TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                    TG2 = (2*X2 - 0.156250000000000000E+01_dp)*0.160000000000000000E+02_dp
                                    CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 162))
                                 ELSE
                                    TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                    TG2 = (2*X2 - 0.168750000000000000E+01_dp)*0.160000000000000000E+02_dp
                                    CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 163))
                                 END IF
                              ELSE
                                 TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 164))
                              END IF
                           ELSE
                              IF (X1 <= 0.312500000000000000E-01_dp) THEN
                                 IF (X2 <= 0.937500000000000000E+00_dp) THEN
                                    IF (X1 <= 0.156250000000000000E-01_dp) THEN
                                       IF (X2 <= 0.906250000000000000E+00_dp) THEN
                                          TG1 = (2*X1 - 0.156250000000000000E-01_dp)*0.640000000000000000E+02_dp
                                          TG2 = (2*X2 - 0.178125000000000000E+01_dp)*0.320000000000000000E+02_dp
                                          CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 165))
                                       ELSE
                                          TG1 = (2*X1 - 0.156250000000000000E-01_dp)*0.640000000000000000E+02_dp
                                          TG2 = (2*X2 - 0.184375000000000000E+01_dp)*0.320000000000000000E+02_dp
                                          CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 166))
                                       END IF
                                    ELSE
                                       TG1 = (2*X1 - 0.468750000000000000E-01_dp)*0.640000000000000000E+02_dp
                                       TG2 = (2*X2 - 0.181250000000000000E+01_dp)*0.160000000000000000E+02_dp
                                       CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 167))
                                    END IF
                                 ELSE
                                    IF (X1 <= 0.156250000000000000E-01_dp) THEN
                                       IF (X2 <= 0.968750000000000000E+00_dp) THEN
                                          IF (X1 <= 0.781250000000000000E-02_dp) THEN
                                             IF (X2 <= 0.953125000000000000E+00_dp) THEN
                                                TG1 = (2*X1 - 0.781250000000000000E-02_dp)*0.128000000000000000E+03_dp
                                                TG2 = (2*X2 - 0.189062500000000000E+01_dp)*0.640000000000000000E+02_dp
                                                CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 168))
                                             ELSE
                                                TG1 = (2*X1 - 0.781250000000000000E-02_dp)*0.128000000000000000E+03_dp
                                                TG2 = (2*X2 - 0.192187500000000000E+01_dp)*0.640000000000000000E+02_dp
                                                CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 169))
                                             END IF
                                          ELSE
                                             TG1 = (2*X1 - 0.234375000000000000E-01_dp)*0.128000000000000000E+03_dp
                                             TG2 = (2*X2 - 0.190625000000000000E+01_dp)*0.320000000000000000E+02_dp
                                             CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 170))
                                          END IF
                                       ELSE
                                          IF (X1 <= 0.781250000000000000E-02_dp) THEN
                                             IF (X2 <= 0.984375000000000000E+00_dp) THEN
                                                TG1 = (2*X1 - 0.781250000000000000E-02_dp)*0.128000000000000000E+03_dp
                                                TG2 = (2*X2 - 0.195312500000000000E+01_dp)*0.640000000000000000E+02_dp
                                                CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 171))
                                             ELSE
                                                TG1 = (2*X1 - 0.781250000000000000E-02_dp)*0.128000000000000000E+03_dp
                                                TG2 = (2*X2 - 0.198437500000000000E+01_dp)*0.640000000000000000E+02_dp
                                                CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 172))
                                             END IF
                                          ELSE
                                             IF (X2 <= 0.984375000000000000E+00_dp) THEN
                                                TG1 = (2*X1 - 0.234375000000000000E-01_dp)*0.128000000000000000E+03_dp
                                                TG2 = (2*X2 - 0.195312500000000000E+01_dp)*0.640000000000000000E+02_dp
                                                CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 173))
                                             ELSE
                                                TG1 = (2*X1 - 0.234375000000000000E-01_dp)*0.128000000000000000E+03_dp
                                                TG2 = (2*X2 - 0.198437500000000000E+01_dp)*0.640000000000000000E+02_dp
                                                CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 174))
                                             END IF
                                          END IF
                                       END IF
                                    ELSE
                                       IF (X2 <= 0.968750000000000000E+00_dp) THEN
                                          TG1 = (2*X1 - 0.468750000000000000E-01_dp)*0.640000000000000000E+02_dp
                                          TG2 = (2*X2 - 0.190625000000000000E+01_dp)*0.320000000000000000E+02_dp
                                          CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 175))
                                       ELSE
                                          IF (X1 <= 0.234375000000000000E-01_dp) THEN
                                             TG1 = (2*X1 - 0.390625000000000000E-01_dp)*0.128000000000000000E+03_dp
                                             TG2 = (2*X2 - 0.196875000000000000E+01_dp)*0.320000000000000000E+02_dp
                                             CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 176))
                                          ELSE
                                             TG1 = (2*X1 - 0.546875000000000000E-01_dp)*0.128000000000000000E+03_dp
                                             TG2 = (2*X2 - 0.196875000000000000E+01_dp)*0.320000000000000000E+02_dp
                                             CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 177))
                                          END IF
                                       END IF
                                    END IF
                                 END IF
                              ELSE
                                 IF (X2 <= 0.937500000000000000E+00_dp) THEN
                                    TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
                                    TG2 = (2*X2 - 0.181250000000000000E+01_dp)*0.160000000000000000E+02_dp
                                    CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 178))
                                 ELSE
                                    IF (X1 <= 0.468750000000000000E-01_dp) THEN
                                       IF (X2 <= 0.968750000000000000E+00_dp) THEN
                                          TG1 = (2*X1 - 0.781250000000000000E-01_dp)*0.640000000000000000E+02_dp
                                          TG2 = (2*X2 - 0.190625000000000000E+01_dp)*0.320000000000000000E+02_dp
                                          CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 179))
                                       ELSE
                                          TG1 = (2*X1 - 0.781250000000000000E-01_dp)*0.640000000000000000E+02_dp
                                          TG2 = (2*X2 - 0.196875000000000000E+01_dp)*0.320000000000000000E+02_dp
                                          CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 180))
                                       END IF
                                    ELSE
                                       TG1 = (2*X1 - 0.109375000000000000E+00_dp)*0.640000000000000000E+02_dp
                                       TG2 = (2*X2 - 0.193750000000000000E+01_dp)*0.160000000000000000E+02_dp
                                       CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 181))
                                    END IF
                                 END IF
                              END IF
                           END IF
                        ELSE
                           IF (X2 <= 0.875000000000000000E+00_dp) THEN
                              TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
                              TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
                              CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 182))
                           ELSE
                              IF (X1 <= 0.937500000000000000E-01_dp) THEN
                                 IF (X2 <= 0.937500000000000000E+00_dp) THEN
                                    TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
                                    TG2 = (2*X2 - 0.181250000000000000E+01_dp)*0.160000000000000000E+02_dp
                                    CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 183))
                                 ELSE
                                    TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
                                    TG2 = (2*X2 - 0.193750000000000000E+01_dp)*0.160000000000000000E+02_dp
                                    CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 184))
                                 END IF
                              ELSE
                                 TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
                                 TG2 = (2*X2 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
                                 CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 185))
                              END IF
                           END IF
                        END IF
                     END IF
                  ELSE
                     IF (X2 <= 0.750000000000000000E+00_dp) THEN
                        TG1 = (2*X1 - 0.375000000000000000E+00_dp)*0.800000000000000000E+01_dp
                        TG2 = (2*X2 - 0.125000000000000000E+01_dp)*0.400000000000000000E+01_dp
                        CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 186))
                     ELSE
                        IF (X1 <= 0.187500000000000000E+00_dp) THEN
                           IF (X2 <= 0.875000000000000000E+00_dp) THEN
                              TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
                              TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
                              CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 187))
                           ELSE
                              TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
                              TG2 = (2*X2 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
                              CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 188))
                           END IF
                        ELSE
                           TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
                           TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
                           CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 189))
                        END IF
                     END IF
                  END IF
               END IF
            ELSE
               IF (X1 <= 0.375000000000000000E+00_dp) THEN
                  IF (X1 <= 0.312500000000000000E+00_dp) THEN
                     IF (X1 <= 0.281250000000000000E+00_dp) THEN
                        TG1 = (2*X1 - 0.531250000000000000E+00_dp)*0.320000000000000000E+02_dp
                        TG2 = (2*X2 - 0.100000000000000000E+01_dp)*0.100000000000000000E+01_dp
                        CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 190))
                     ELSE
                        TG1 = (2*X1 - 0.593750000000000000E+00_dp)*0.320000000000000000E+02_dp
                        TG2 = (2*X2 - 0.100000000000000000E+01_dp)*0.100000000000000000E+01_dp
                        CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 191))
                     END IF
                  ELSE
                     IF (X2 <= 0.500000000000000000E+00_dp) THEN
                        TG1 = (2*X1 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
                        TG2 = (2*X2 - 0.500000000000000000E+00_dp)*0.200000000000000000E+01_dp
                        CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 192))
                     ELSE
                        TG1 = (2*X1 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
                        TG2 = (2*X2 - 0.150000000000000000E+01_dp)*0.200000000000000000E+01_dp
                        CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 193))
                     END IF
                  END IF
               ELSE
                  IF (X1 <= 0.437500000000000000E+00_dp) THEN
                     IF (X2 <= 0.500000000000000000E+00_dp) THEN
                        TG1 = (2*X1 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
                        TG2 = (2*X2 - 0.500000000000000000E+00_dp)*0.200000000000000000E+01_dp
                        CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 194))
                     ELSE
                        IF (X1 <= 0.406250000000000000E+00_dp) THEN
                           TG1 = (2*X1 - 0.781250000000000000E+00_dp)*0.320000000000000000E+02_dp
                           TG2 = (2*X2 - 0.150000000000000000E+01_dp)*0.200000000000000000E+01_dp
                           CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 195))
                        ELSE
                           TG1 = (2*X1 - 0.843750000000000000E+00_dp)*0.320000000000000000E+02_dp
                           TG2 = (2*X2 - 0.150000000000000000E+01_dp)*0.200000000000000000E+01_dp
                           CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 196))
                        END IF
                     END IF
                  ELSE
                     IF (X2 <= 0.500000000000000000E+00_dp) THEN
                        TG1 = (2*X1 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
                        TG2 = (2*X2 - 0.500000000000000000E+00_dp)*0.200000000000000000E+01_dp
                        CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 197))
                     ELSE
                        TG1 = (2*X1 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
                        TG2 = (2*X2 - 0.150000000000000000E+01_dp)*0.200000000000000000E+01_dp
                        CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 198))
                     END IF
                  END IF
               END IF
            END IF
         ELSE
            IF (X1 <= 0.750000000000000000E+00_dp) THEN
               IF (X1 <= 0.625000000000000000E+00_dp) THEN
                  IF (X2 <= 0.500000000000000000E+00_dp) THEN
                     IF (X1 <= 0.562500000000000000E+00_dp) THEN
                        TG1 = (2*X1 - 0.106250000000000000E+01_dp)*0.160000000000000000E+02_dp
                        TG2 = (2*X2 - 0.500000000000000000E+00_dp)*0.200000000000000000E+01_dp
                        CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 199))
                     ELSE
                        TG1 = (2*X1 - 0.118750000000000000E+01_dp)*0.160000000000000000E+02_dp
                        TG2 = (2*X2 - 0.500000000000000000E+00_dp)*0.200000000000000000E+01_dp
                        CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 200))
                     END IF
                  ELSE
                     IF (X1 <= 0.562500000000000000E+00_dp) THEN
                        TG1 = (2*X1 - 0.106250000000000000E+01_dp)*0.160000000000000000E+02_dp
                        TG2 = (2*X2 - 0.150000000000000000E+01_dp)*0.200000000000000000E+01_dp
                        CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 201))
                     ELSE
                        TG1 = (2*X1 - 0.118750000000000000E+01_dp)*0.160000000000000000E+02_dp
                        TG2 = (2*X2 - 0.150000000000000000E+01_dp)*0.200000000000000000E+01_dp
                        CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 202))
                     END IF
                  END IF
               ELSE
                  IF (X2 <= 0.500000000000000000E+00_dp) THEN
                     IF (X1 <= 0.687500000000000000E+00_dp) THEN
                        TG1 = (2*X1 - 0.131250000000000000E+01_dp)*0.160000000000000000E+02_dp
                        TG2 = (2*X2 - 0.500000000000000000E+00_dp)*0.200000000000000000E+01_dp
                        CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 203))
                     ELSE
                        TG1 = (2*X1 - 0.143750000000000000E+01_dp)*0.160000000000000000E+02_dp
                        TG2 = (2*X2 - 0.500000000000000000E+00_dp)*0.200000000000000000E+01_dp
                        CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 204))
                     END IF
                  ELSE
                     TG1 = (2*X1 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
                     TG2 = (2*X2 - 0.150000000000000000E+01_dp)*0.200000000000000000E+01_dp
                     CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 205))
                  END IF
               END IF
            ELSE
               IF (X1 <= 0.875000000000000000E+00_dp) THEN
                  TG1 = (2*X1 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
                  TG2 = (2*X2 - 0.100000000000000000E+01_dp)*0.100000000000000000E+01_dp
                  CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 206))
               ELSE
                  TG1 = (2*X1 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
                  TG2 = (2*X2 - 0.100000000000000000E+01_dp)*0.100000000000000000E+01_dp
                  CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 207))
               END IF
            END IF
         END IF
      END IF
   END SUBROUTINE t_c_g0_n

! **************************************************************************************************
!> \brief ...
!> \param Nder the number of derivatives that will actually be used
!> \param iunit contains the data file to initialize the table
!> \param mepos ...
!> \param group ...
! **************************************************************************************************
   SUBROUTINE init(Nder, iunit, mepos, group)
      INTEGER, INTENT(IN)                                :: Nder, iunit, mepos

      CLASS(mp_comm_type), INTENT(IN)                     :: group

      INTEGER                                            :: I
      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: chunk

      patches = 207
      IF (Nder > nderiv_max) CPABORT("T_C_G0 init failed")
      nderiv_init = Nder
      IF (ALLOCATED(C0)) DEALLOCATE (C0)
      ! round up to a multiple of 32 to give some generous alignment for each C0
      ALLOCATE (C0(32*((31 + (Nder + 1)*(degree + 1)*(degree + 2)/2)/32), patches))
      ! init mpi'ed buffers to silence warnings under valgrind
      C0 = 1.0E99_dp
      IF (mepos == 0) THEN
         ALLOCATE (chunk((nderiv_max + 1)*(degree + 1)*(degree + 2)/2))
         DO I = 1, patches
            READ (iunit, *) chunk
            C0(1:(Nder + 1)*(degree + 1)*(degree + 2)/2, I) = chunk(1:(Nder + 1)*(degree + 1)*(degree + 2)/2)
         END DO
         DEALLOCATE (chunk)
      END IF
      CALL group%bcast(C0, 0)

   END SUBROUTINE init

! **************************************************************************************************
!> \brief ...
! **************************************************************************************************
   SUBROUTINE free_C0()
      IF (ALLOCATED(C0)) DEALLOCATE (C0)
      nderiv_init = -1
   END SUBROUTINE free_C0

! **************************************************************************************************
!> \brief ...
!> \param RES ...
!> \param NDERIV ...
!> \param TG1 ...
!> \param TG2 ...
!> \param C0 ...
! **************************************************************************************************
   SUBROUTINE PD2VAL(RES, NDERIV, TG1, TG2, C0)
      REAL(KIND=dp), INTENT(OUT)                         :: res(*)
      INTEGER, INTENT(IN)                                :: NDERIV
      REAL(KIND=dp), INTENT(IN)                          :: TG1, TG2, C0(105, *)

      REAL(KIND=dp), PARAMETER :: SQRT2 = 1.4142135623730950488016887242096980785696718753_dp

      INTEGER                                            :: K
      REAL(KIND=dp)                                      :: T1(0:13), T2(0:13)

      T1(0) = 1.0_dp
      T2(0) = 1.0_dp
      T1(1) = SQRT2*TG1
      T2(1) = SQRT2*TG2
      T1(2) = 2*TG1*T1(1) - SQRT2
      T2(2) = 2*TG2*T2(1) - SQRT2
      T1(3) = 2*TG1*T1(2) - T1(1)
      T2(3) = 2*TG2*T2(2) - T2(1)
      T1(4) = 2*TG1*T1(3) - T1(2)
      T2(4) = 2*TG2*T2(3) - T2(2)
      T1(5) = 2*TG1*T1(4) - T1(3)
      T2(5) = 2*TG2*T2(4) - T2(3)
      T1(6) = 2*TG1*T1(5) - T1(4)
      T2(6) = 2*TG2*T2(5) - T2(4)
      T1(7) = 2*TG1*T1(6) - T1(5)
      T2(7) = 2*TG2*T2(6) - T2(5)
      T1(8) = 2*TG1*T1(7) - T1(6)
      T2(8) = 2*TG2*T2(7) - T2(6)
      T1(9) = 2*TG1*T1(8) - T1(7)
      T2(9) = 2*TG2*T2(8) - T2(7)
      T1(10) = 2*TG1*T1(9) - T1(8)
      T2(10) = 2*TG2*T2(9) - T2(8)
      T1(11) = 2*TG1*T1(10) - T1(9)
      T2(11) = 2*TG2*T2(10) - T2(9)
      T1(12) = 2*TG1*T1(11) - T1(10)
      T2(12) = 2*TG2*T2(11) - T2(10)
      T1(13) = 2*TG1*T1(12) - T1(11)
      T2(13) = 2*TG2*T2(12) - T2(11)
      DO K = 1, NDERIV + 1
         RES(K) = 0.0_dp
         RES(K) = RES(K) + DOT_PRODUCT(T1(0:13), C0(1:14, K))*T2(0)
         RES(K) = RES(K) + DOT_PRODUCT(T1(0:12), C0(15:27, K))*T2(1)
         RES(K) = RES(K) + DOT_PRODUCT(T1(0:11), C0(28:39, K))*T2(2)
         RES(K) = RES(K) + DOT_PRODUCT(T1(0:10), C0(40:50, K))*T2(3)
         RES(K) = RES(K) + DOT_PRODUCT(T1(0:9), C0(51:60, K))*T2(4)
         RES(K) = RES(K) + DOT_PRODUCT(T1(0:8), C0(61:69, K))*T2(5)
         RES(K) = RES(K) + DOT_PRODUCT(T1(0:7), C0(70:77, K))*T2(6)
         RES(K) = RES(K) + DOT_PRODUCT(T1(0:6), C0(78:84, K))*T2(7)
         RES(K) = RES(K) + DOT_PRODUCT(T1(0:5), C0(85:90, K))*T2(8)
         RES(K) = RES(K) + DOT_PRODUCT(T1(0:4), C0(91:95, K))*T2(9)
         RES(K) = RES(K) + DOT_PRODUCT(T1(0:3), C0(96:99, K))*T2(10)
         RES(K) = RES(K) + DOT_PRODUCT(T1(0:2), C0(100:102, K))*T2(11)
         RES(K) = RES(K) + DOT_PRODUCT(T1(0:1), C0(103:104, K))*T2(12)
         RES(K) = RES(K) + DOT_PRODUCT(T1(0:0), C0(105:105, K))*T2(13)
      END DO
   END SUBROUTINE PD2VAL

! **************************************************************************************************
!> \brief Returns the value of nderiv_init so that one can check if opening the potential file is
!>        worhtwhile
!> \return ...
!> \author A. Bussy, 05.2019
! **************************************************************************************************
   FUNCTION get_lmax_init() RESULT(lmax_init)

      INTEGER                                            :: lmax_init

      lmax_init = nderiv_init

   END FUNCTION get_lmax_init

END MODULE t_c_g0
